Age group DNA methylation differences in lemon sharks (Negaprion brevirostris): Implications for future age estimation tools

Abstract Age information is often non‐existent for most shark populations due to a lack of measurable physiological and morphological traits that can be used to estimate age. Recently, epigenetic clocks have been found to accurately estimate age for mammals, birds, and fish. However, since these clocks rely, among other things, on the availability of reference genomes, their application is hampered in non‐traditional model organisms lacking such molecular resources. The technique known as Methyl‐Sensitive Amplified Polymorphism (MSAP) has emerged as a valid alternative for studying DNA methylation biomarkers when reference genome information is missing, and large numbers of samples need to be processed. Accordingly, the MSAP technique was used in the present study to characterize global DNA methylation patterns in lemon sharks from three different age groups (juveniles, subadults, and adults). The obtained results reveal that, while MSAP analyses lack enough resolution as a standalone approach to infer age in these organisms, the global DNA methylation patterns observed using this technique displayed significant differences between age groups. Overall, these results confer that DNA methylation does change with age in sharks like what has been seen for other vertebrates and that MSAP could be useful as part of an epigenetics pipeline to infer the broad range of ages found in large samples sizes.


| INTRODUC TI ON
Overfishing constitutes the most important threat to elasmobranchs including sharks, skates, and rays (Dulvy et al., 2014), making fisheries management critical for their conservation. Consequently, the characterization of demographic information is of utmost importance to model the sustainability of populations under fishing pressure. Age information is a critical parameter of these models, which is unfortunately not available for most species due (among other reasons) to the lack of physiological and morphological traits that can be accurately used as proxies of chronological age in elasmobranchs. Particularly in the case of sharks, body length has been largely used to infer age and has often been cross-validated with vertebral growth bands (Carlson & Goldman, 2007). Yet, recent studies have reported that the age of older sharks has been vastly underestimated given that vertebral bands do not necessarily form on a yearly basis throughout a shark's lifetime (Natanson et al., 2018;Natanson & Skomal, 2015). Not less importantly, vertebral band counts are highly invasive, require sacrificing animals, and do not work for all species (Francis et al., 2007;Huveneers et al., 2013;Natanson et al., 2018).
The emergence of methodological approaches linking molecular variation to chronological age has constituted a breakthrough, particularly in those cases where age estimation has proven difficult and/or based on invasive methods. The characterization of distinctive molecular traits is possible using small tissue samples that can be more easily obtained, in most cases without even the need to capture/release the animal. One of the first among these methods was based on the quantification of the length of telomeres, based on their progressive shortening during successive cell cycle divisions.
While telomere-based aging methods proved to be a great proxy for biological age (i.e., how old an animal looks for its age), the vast differences in the starting length observed in telomeres led numerous studies to question its applicability for determining chronological age estimation (i.e., actual age of the organism) (Nehmens et al., 2021;Remot et al., 2020;Vaiserman & Krasnienkov, 2020).
Epigenetics, defined as "the study of molecules and mechanisms that can perpetuate alternative gene activity states in the context of the same DNA sequence" (Cavalli & Heard, 2019), has provided a basis to start solving that problem by establishing direct links between chronological age and particular epigenetic modifications including DNA methylation, histone variants, histone modifications, and small RNAs (López-Otín et al., 2013), providing not only insights into how organisms age but also facilitating the estimation of age using DNA methylation (the addition of a methyl group most commonly to a cytosine in the DNA sequence) as a proxy (Horvath, 2013).
Chronological epigenetic clocks specifically correlate percent DNA methylation at CpG sites (a Cytosine followed by a Guanine in the same DNA strand) to age through multiple regression analyses producing highly accurate estimations (correlation coefficient, r 2 , ranging from 0.70 to 0.99). These changes in DNA methylation are linked to age-related mechanisms and developmental milestones that often result in changes in gene expression (Gladyshev, 2016;Hannum et al., 2013;López-Otín et al., 2013). Although largely focused on mammalian species (Barratclough et al., 2021;Beal et al., 2019;Bors et al., 2021;Fei et al., 2021;Horvath & Raj, 2018;Polanowski et al., 2014;Robeck et al., 2021;Stubbs et al., 2017;Tanabe et al., 2020), accurate clocks have been developed for other taxa including birds and fishes (Anastasiadi & Piferrer, 2020 Chronological epigenetic clocks get their accuracy from the characterization of DNA methylation at particular DNA positions using high-resolution approaches such as microarrays, reduced representation bisulfite sequencing, pyrosequencing, and qPCR (Morselli et al., 2021;Trigg et al., 2021). However, such techniques are hampered by two major limitations: the requirement of a sequenced reference genome for the species studied and the potentially high cost of conducting these analyses in large numbers of samples to be ecologically meaningful in field studies, particularly for the larger sample sizes necessary for ecologically relevant population and stock assessments. Thus, the application of less demanding and inexpensive coarse-grained methods, still providing useful aging information at the level of age intervals, represents a compromise helping the characterization of aging structures in natural populations within realistic ecological settings. The characterization of Methyl-Sensitive Amplified Polymorphism (MSAP) is one such method, based on the use of dual methylation-sensitive restriction enzymes to unveil broad DNA methylation patterns found at target cut site 5'-CCGG-3′ (invented by Reyna-López et al., 1997). This method has been successfully used to examine the links between environmental change and genome-wide DNA methylation patterns involved in phenotypic responses (Beal et al., 2021;Huang et al., 2017;Morán et al., 2013;Pierron et al., 2014;Rodríguez-Casariego et al., 2020;Suarez-Ulloa et al., 2019;Zhao et al., 2015).
Since it is well known that gene expression patterns change during development and aging, the observation of global amounts of DNA methylation changing over an individual's lifetime comes as no surprise (Heyn et al., 2012). However, it is not yet known if such changes are large enough to be efficiently detected using MSAP as a standalone epigenetic aging method. The present study addresses that question by using MSAP in lemon shark (Negaprion brevirostris) samples belonging to three age groups (juvenile, subadult, and adult), available through the unique long-term sampling archive from Bimini, Bahamas. This archive of samples is one of the most comprehensive datasets for sharks and is the best known age dataset to date for any shark species. The MSAP method used in the present study was previously adjusted for use in lemon sharks from populations monitored in our own research (Beal et al., 2021). The obtained results are consistent with the applicability of MSAP analyses as part of a larger epigenetic pipeline to determine chronological age, particular at earlier stages, to define age groups in large numbers of samples which can be subsequently more finely characterized using single nucleotide resolution methods.

| Samples and DNA extraction
Lemon sharks, as studied in Bimini (Bahamas) by the Bimini Shark Lab for more than 30 years, exhibit a life history of biennial parturition after reaching sexual maturity around age 12 (Brown & Gruber, 1988;Feldheim et al., 2004) as well as natal philopatry where they return back to the particular nursery site where they were born (Feldheim et al., 2014). Although age of reproductive senescence is unknown, pregnant females in the low to mid-30s have been recorded in the Bimini area and the maximum lifespan is estimated at 37 years (Brooks et al., 2016).
Epigenetic age estimation has been successful across multiple tissue types (Horvath, 2013). Therefore, we focused on fin clips as this is the type of tissue sample most often collected during shark population monitoring efforts and therefore most relevant for managers and conservation practitioners. All lemon shark samples analyzed in this study consisted of fin clips housed in an archived collection at the Field Museum in Chicago with collection dates ranging Subadults (6-11 years) were between 70-175 cm PCL for males and 70-185 cm for females, and adults (12+ years) were greater than 175 cm PCL for males and greater than 185 cm for females.
The samples used in this study were processed in two batches.
The first batch (Batch 1, LS1-10) was analyzed as part of a proof-ofconcept pilot study, where DNA extraction was performed at the Field Museum in Chicago using a salting-out protocol (Sunnucks & Hales, 1996). These samples were preserved in ice and shipped to Florida International University (FIU) for MSAP analyses. Batch 2 samples were shipped as fin clip samples in DMSO, and DNA was subsequently extracted following the same salting-out protocol used for Batch 1 samples (Beal et al., 2021). Batch 1 samples (n = 10) consisted of four juveniles, three subadults, and three adults. Batch 2 samples (n = 19) consisted of six juveniles, nine subadults, and four adults.

| DNA methylation pattern analysis
DNA methylation patterns were determined using MSAP analyses. This method employs dual methylation-sensitive enzymes that cut differentially at the same target site 5'-CCGG-3′ (Reyna-López et al., 1997) creating dual fragment profiles identifying four types of DNA methylation states present (full or hypermethylation, internal C methylation, hemi-methylation, or no methylation).
MSAP analyses were conducted independently for each of the two batches of samples as detailed in Beal et al. (2021). Briefly, a parallel restriction digestion/ligation (each containing one of the methylsensitive restriction enzymes MSPI or HPAII and each containing the non-methyl sensitive rare cutting restriction enzyme EcoRI) was performed using 200 ng of DNA, followed by a pre-selective PCR targeting the ligated ends with an additional base pair, and ending with a final selective PCR with an additional two base pairs added to the primer along with a fluorescent dye (6-FAM) used to detect fragment length (Table 2). A total of four selective primers in two combinations were used for both batches. Fragment analysis was conducted by the FIU DNA Core on an Applied Biosystems® 3130XL Genetic Analyzer. Raw fragment data were then scored as presence (+)/ absence (−) for all fragments for each individual for each parallel reaction (MSPI and HPAII) to decipher the four types of DNA methylation status present. Presence (+)/absence (−) (MSPI/HPAII) pattern for each fragment gave methylation type (+/+ = unmethylated, +/− = hemi-methylated, −/+ = internal C methylation, −/− = full methylation). Batch 1 samples were subjected to the MSAP protocol and fragment analysis at a separate and earlier time than Batch 2 samples. Once raw data were available for all samples, the final dataset was scored using all samples to allow for loci to be compared.

| Data analysis
The msap package (version 1.1.8) in R was used to categorize each locus as either unmethylated, hemi-methylated, methylated at an internal cytosine, or fully methylated, for each sample (Pérez-Figueroa, 2013). Although a lack of bands in both reactions could indicate full methylation or the lack of a genetic target, this case was considered as full methylation, consistent with other DNA methylation studies working with species where there is low genetic structure (Feldheim et al., 2001;Rodríguez-Casariego et al., 2020), as has been shown in lemon sharks (Feldheim et al., 2001). Methylationsusceptible loci (MSL) were identified as loci with a methylated state in at least 5% of the samples (Pérez-Figueroa, 2013), and polymorphic loci were identified as those that displayed both methylation and a lack of methylation in at least two instances across samples (Herrera & Bazaga, 2010;Pérez-Figueroa, 2013) The optimal number of principal components (PCs) to be retained (4 PCs) was determined through the xval-Dapc function as well as an assessment of a-scores (optim.a.score function; Jombart et al., 2010).
Both discriminant functions were retained. To assess the separation between age groups, the probability of membership of each sample to each age group was evaluated. A distance-based redundancy analysis (dbRDA; Legendre & Anderson, 1999;Oksanen et al., 2019) was performed with batch 2 samples to characterize the relationship between all the metadata available for the samples (age class, sex, year, season) and DNA methylation variation. The equation was as follows: A permutation test was performed on the dbRDA model (anova. cca function) to assess the significance of these variables upon the model (Legendre et al., 2011).

| Combined batch sample analysis
A preliminary MSAP analysis using Batch 1 samples (n = 10) revealed significant differences between age groups (Analysis of molecular variance (AMOVA), p < .05, Figures S1 and S2). In order to validate these findings, additional samples were incorporated into the analyses (Batch 1 + Batch 2 samples). The combined batch analysis encompassed a total of n = 29 samples (10 juveniles, 12 subadults, and seven adults) used in R analyses, specifically the msap package (Table 1). Two selective primer pairs (Table 2) were used which resulted in a total of ( ) ∼ + + + . and Batch 2 (Figure 1). The potential effect of differences in year of sampling and season was investigated, but neither appeared to have a significant influence in the clustering of data points ( Figures S3 and   S4). The only obvious difference appeared to be the fact that the samples were processed in two separate batches that were not processed by the same lab for all steps. Therefore, it was decided to investigate age group clustering within each batch effort to eliminate interference from the different sample processing and storage variables. Despite this apparent batch effect, age group differences were readily obvious and appeared robust despite possible differences between batches.

| Independent batch sample analysis
The analysis of Batch 1 samples (n = 10) consisted of four juveniles, three subadults, and three adults. A total of 308 loci were identified as MSL (45% of them polymorphic) and 56 to as NML (0% polymorphic) using the msap package ( Figures S5 and S6). AMOVA analyses revealed a significant difference between age groups for MSL (Phi_ST = 0.111, df = 9, p = .049). Accordingly, juveniles had less unmethylated loci compared to subadult and adult groups (Figure 2).
The NMDS plot showed that juveniles cluster more closely together than individuals in subadult or adult groups (Figure 3). Similarly, PERMANOVA analyses found a significant difference between age groups (p = .05) and a marginally significant difference in dispersion between groups (p = .09). Such age group differences were similar to those observed in the combined batch sample analysis, with samples

TA B L E 2 Adaptors and Primers usedfor MSAP analyses. Only selective primers were fluorescently labeled (FAM) for fragment analysis
Step/combo Oligo name Sequence (5′ to 3′)  F I G U R E 3 NMDS plot for batch 1. Distance between any two dots depicts the difference/similarity between two given samples F I G U R E 4 Percent DNA methylation per age Group for Batch 2. Plot (a) depicts % loci as methylated or non-methylated. Plot (b) depicts % loci for each of the four types of methylation status determined by MSAP analysis (none-non-methylated, hemi-hemi-methylated, Internal C-internal cytosine methylated, full-full/hypermethylated target site).

Restriction ligation EcoRI Adaptor Fwd CTCGTAGACTGCGTACC
from juvenile sharks clustering independently from subadults and adults ( Figure 5). PERMANOVA analysis identified significant differences between age groups and dispersion between age groups (p = .001 and p = .007, respectively). To further analyze differences between age groups, a DAPC analysis was conducted on Batch 2 samples (feasible in this particular case based on its larger sample size) with discriminant function 1 accounting for the majority (~50%) of differentiation observed between juveniles and other age groups ( Figure 6). Lemon shark age groups (i.e., juvenile, subadult, and adult) showed differential DNA methylation patterns ( Figures 5 and 6), with the juvenile age group showing the most differentiated pattern from the other groups and with 100% accurate membership predictability compared with individuals from the other two age groups (Figure 7).
These results are consistent with a progressive change in DNA methylation pattern from juvenile to adult individuals (Figure 6b) F I G U R E 5 NMDS plot for batch 2. Distance between any two dots depicts the difference/similarity between two given samples F I G U R E 6 DAPC and density plots for discriminant function 1 and 2 for batch 2. Plot (a), shows spatial and temporal distribution of data points for the three age groups, x-axis is DF1 and y-axis is DF2. Plot (b), density plot for DF1 showing the most differentiation of juveniles from other age groups. Plot (c), density plot for DF2. In both B and C plots, we see that subadults have the most variation across both DF1 and DF2.

| DISCUSS ION
This study represents the first effort to investigate the relation- Beyond the observed batch effect, the differences by age group were visible for combined Batch (1 + 2) analysis (n = 29; Figure 1) as well as when each batch (Batch 1, n = 10, Figure 3; and Batch 2, n = 19, Figure 5) was analyzed individually. The consistency of the age group pattern regardless of the batch analyzed further supports a strong link between age and DNA methylation patterns across age groups. To evaluate the observed age group differences, further analyses were conducted using batch 2 samples only in order to remove any interference from the batch effect (Figures 6 and 7) based on the larger sample size of this batch (n = 19).

| The importance of age clustering for clock development
The findings of this study provide the first evidence for the role of epigenetics during aging in elasmobranchs including sharks, a mysterious group known for the long life spans in some of its members.
Previous studies on the whale shark reported that global DNA methylation patterns were similar to other vertebrate species (i.e., high levels of methylation at most CpG sites except for those at promoter regions), leading to the conclusion that this epigenetic modification plays a similar role in sharks as that characterized in other vertebrates (Peat et al., 2017). The present study further supports that conclusion, based on the similarity between the observed age-related changes in DNA methylation patterns in lemon sharks and other vertebrates (Anastasiadi & Piferrer, 2020;Beal et al., 2019;Bors et al., 2021;De Paoli-Iseppi et al., 2019;Fei et al., 2021;Polanowski et al., 2014;Stubbs et al., 2017). The largest age-related change in DNA methylation patterns was found between juveniles and adults seen in the DAPC analysis (Figure 6a), with approximately 50% of the differences between groups explained by discriminant function 1 (x-axis Figure 6a,b). Interestingly, subadult DNA methylation patterns had a large amount of variation that spanned between the juvenile and adult groups (Figure 6b). These results seem to be in line with reports describing gradual overall changes in DNA methylation throughout the genome correlated with chronological aging (Bollati et al., 2009;Ciccarone et al., 2018;Day et al., 2013;Horvath, 2013;Maegawa et al., 2010). Along the continuum of age, hypomethylation has been observed across the genome due to drift and errors in maintenance of DNA methylation on the genome (Berdyshev et al., 1967;Horvath, 2013;Martin-Herranz et al., 2019;Shimoda et al., 2014;Wilson et al., 1987). However, this is not the case for all sites as revealed by the presence of some sites across the genome displaying age-driven hypermethylation (Rakyan et al., 2010).

F I G U R E 7
Probability of membership to a given age group. Samples are listed on the x-axis as sample ID_Letter representative of age group as follows: J-juveniles, S-subadult, and a-adult. Juveniles showed a 100% probability of being assigned to their correct group. Subadults had the most variation (50%-95%) in probability of correct assignment, whereas all adult samples had 60% or greater probability.
Some of these changes are possibly correlated with changes in gene expression during development (Calvanese et al., 2009;López-Otín et al., 2013;Pal & Tyler, 2016;Unnikrishnan et al., 2019). This bidirectional change along with the specificity of individual CpG sites correlating to age could explain why subadults, as observed in this study, have a wider variation in DNA methylation points that seems to span between juveniles and adults as they transition between these two age groups. It is also possible that this same bi-directional change is one reason that the differentiation between subadults and adults is less obvious and may require more higher resolution methods to elucidate age for these groups.
Another interesting pattern observed in this work is the larger amount of variation in DNA methylation patterns for older shark individuals. Many of these DNA methylation positions most likely are driven by the onset of epigenetic modifications linked to changes in gene expression related to aging; however, many may very well be attributed to the higher incidence of epigenetic drift later in life associated with heterogeneous environments (Jung et al., 2017), contributing to the variation observed, especially since these samples come from varying sample years. This effect is best illustrated by human twin studies revealing a variation in epigenomes driven by environmental heterogeneity (Bell & Spector, 2011;Fraga et al., 2005) but has also been captured by epigenetic clock studies (Beal et al., 2019;Robeck et al., 2021;Stubbs et al., 2017;Thompson et al., 2017).
Although the results reported in this study support a place for MSAP in aging studies, the prevalence of such random and environmentally driven DNA methylation changes during aging prevents this technique from producing accurate estimations of chronological age. On the other hand, MSAP proved useful for making general assignments of age group. Furthermore, these findings suggest that the epigenetic aging pattern observed in mammals also extends to sharks, validating the widespread role of DNA methylation in such process.

| The manifold implications of the age of the shark
Future development of epigenetic clocks in sharks would benefit shark conservation similar to the benefits that have been found in being able to characterize the direct link between age, environmental stimulus, and age-related diseases developed in humans (Bektas et al., 2018;Cavalli & Heard, 2019;Dugué et al., 2018;Zhang et al., 2017). Better understanding of this connection in sharks would allow for environmental threats impacting populations to be identified, as well as the specific age groups most vulnerable to the threats, something that is little known for most shark species. Additionally, with little to no age data existing for most shark populations, little demographic and life history information is available, and the status of populations is hard to assess (Carlson & Goldman, 2007;Goldman & Cailliet, 2004). The current methods used to age sharks (band central staining and bomb radiocarbon aging) are not applicable to all species, require sacrificing individuals, and are fraught with high variation (Campana, 2001;Francis et al., 2007;Huveneers et al., 2013;Natanson et al., 2018).
Furthermore, being able to estimate shark age from tissue would greatly enhance the conservation work being done in the fin trade which has been identified as one of the top threats to endangered shark species (Walls & Dulvy, 2020). Current molecular methods only allow for the identification of species and possibly the origin population (Cardeñosa et al., 2020(Cardeñosa et al., , 2021Clarke et al., 2004Clarke et al., , 2006Fields et al., 2015Fields et al., , 2018 as well as investigations into the quantities of metals and other pollutants (Barcia et al., 2020;Nalluri et al., 2014) these sharks are exposed to but do not allow for the estimation of age.
Epigenetic alteration is one of the hallmarks of aging, and the study of these alterations could lead to a better understanding and ability to combat age-related diseases (López-Otín et al., 2013).
As such, studying the specific markers that correlate with age for sharks may prove to be beneficial to understanding and learning how to treat age-related diseases in humans. Many sharks live long lives Natanson & Skomal, 2015;Nielsen et al., 2016), and one means to their longevity is thought to be that sharks are less susceptible to age-related diseases. For example, sharks may be resistant to the age-related cancers due to potential tumor-suppressing properties found in their cartilage (although research remains inconclusive) (Oliveira et al., 2019;Walsh et al., 2013); however, there remains a lack of systematic surveys and studies to confirm the actual incidence of cancer in sharks in the first place (Ostrander et al., 2004). Findings from the white shark genome do provide evidence that if resistance to cancer does exist for sharks, it may be due to positive selection of genome stabilizing genes which would have a direct effect upon tumor suppression since many tumors form due to breaks and damage to the DNA (Marra et al., 2019). Further study into epigenetic clocks may provide more insights into the longevity many sharks have and could provide insights into human health and aging.

| CON CLUS IONS
This study contributes to fill two major gaps in the study of shark aging and epigenetics: sharks do exhibit DNA methylation patterns that change with age, and the MSAP method is efficient in identifying such changes. These findings are important because they support the development of epigenetic aging tools in this group of elasmobranchs, fostering the development of additional resources including the generation of full reference genomes. In addition, this work underscores the role of MSAP as a simple, inexpensive, and widely applicable technique to perform the first round of analyses within more elaborated epigenetic pipelines able to systematically estimate chronological age in samples from diverse origins and age groups. Overall, this study provides a first step toward the development of such pipelines for the shark research community with critical implications in conservation and management (Marra et al., 2019).